Rheology of a Supercooled Polymer Melt 
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Molecular dynamics simulations are performed for a polymer melt composed of short chains 
in quiescent and sheared conditions. The stress relaxation function G(t) exhibits a stretched ex- 
ponential form in a relatively early stage and ultimately follows the Rouse function in quiescent 
supercooled state. Transient stress evolution after application of shear obeys the linear growth 
J. dt'G(t') for strain less than 0.1 and then saturates into a non- Newtonian viscosity. In steady 
states, strong shear-thinning and elongation of chains into ellipsoidal shapes are found at extremely 
small shear. A glassy component of the stress is much enhanced in these examples. 

PACS numbers: 83.10.Nn, 83.20.Jp, 83.50.By, 64.70.Pf 



Stress and dielectric relaxations of glassy polymer 
melts occur from microscopic to macroscopic time scales 
in very complicated manners p],@]. Experiments have 
shown that the stress relaxation function G(t) exhibits a 
glassy stretched exponential decay, a glass-rubber tran- 
sition, a rubbery plateau, and a terminal decay, in this 
order over many decades of time. Such hierarchical re- 
laxation behavior arises from rearrangements of jammed 
atomic configurations and subsequent evolution of chain 
conformations described by the Rouse or reptation dy- 
namics Pfl]. The stress-optical relation between bire- 
fringence and stress has also been reported to be violated 
as the temperature T is approached the glass transition 
temperature T g g-Q], obviously owing to enhancement 
of a glassy part of the stress. 

Recent simulations on glassy polymer melts have 
mainly treated self-motions of particles in quiescent 
states [p|-pO[. However, not enough theoretical efforts 
have been made on the rheological properties of glassy 
polymers. Hence, we will first study linear rheology 
of a model short chain system via very long molecular 
dynamics simulations. Then, we will demonstrate that 
chains are very easily elongated at extremely small shear 
rate j on the order of the inverse Rouse time. Marked 
shear-thinning behavior then takes place for larger shear 
rates, indicating a decrease of the monomeric friction 
among different chains. On the other hand, in super- 
cooled simple fluid mixtures pd| , the shear-dependent 
structural rearrangement time 11,(7) depends on shear as 
T b(7) _1 = T &(0) * + Abj, where A\, is of order 1 and 
Tft(O) is on the order of the so-called a relaxation time T a 
obtained from the incoherent van Hove time correlation 
function. The steady state viscosity is expressed in a very 
simple form, 77(7) = A 7] Tb{ A f)-\-T)Bi for any T and 7, where 
A v and t\b are constants. We have also found that the 
cage breakage occurs collectively in the form of clusters 
characterized by a correlation length £ |12J, where the 
dynamic scaling 77, (7) ~ £ 2 holds in three dimensions. 

In our model all the bead particles interact with a 
Lennard-Jones potential of the form [G]||]l(J, Ulj(t) = 



4e[(a/r) 12 — (a/r) 6 ] + e. It is cut off at the minimum 
distance 2 1 ' 6 ct, so we use its repulsive part only to pre- 
vent spatial overlap of particles. Consecutive beads on 
each chain are connected by an anharmonic spring of the 
form, U F (r) = -ifc c i?gln[l - {r/R Q ) 2 } with k c = 30e/V 
and R = 1.5a, so the bond length cannot exceed R . 
In a cubic box with length L — 10a under the periodic 
boundary condition, we put M = 100 chains composed of 
N = 10 beads. The number density is fixed at a very high 
value of n = NM/V = 1/a 3 , which results in severely 
jammed configurations at low T. We will measure space 
and time in units of a and tq = (ma 2 /e) 1 ' 2 with m being 
the mass of a bead. The temperature T will be measured 
in units of e/fc^- Simulations were performed in normal 
(T = 1.0) and supercooled (T = 0.4 and 0.2) states with 



and without shear. The bond lengths bj 



\Rj ~ -ftj+i 



between consecutive beads on each chain exhibit only 
small deviations on the order of a few % from b = 0.96, 
which gives the minimum of Ul,i(t) + Up{r), for any T 
and 7 in our study. Thus the bond lengths are nearly 
fixed in our model system. 

We took data after long equilibration periods (10 6 at 
T = 0.2) to suppress aging (slow equilibration) effects 
in various quantities such as the pressure or the density 
time correlation functions. At zero shear we imposed the 
micro-canonical condition with the time step At = 0.005. 
In order to obtain accurate linear viscoelastic behavior, 
very long simulations of order 10 2 tr were performed, 
where tr is the primary Rouse relaxation time. In the 
literature |4|JgjJlO||, simulation times have been typically 
up to tr in supercooled states. In the presence of shear 
we set At = 0.0025 and kept the temperature at a con- 
stant using the Gaussian constraint thermostat to elim- 
inate viscous heating. After a long equilibration time in 
a quiescent state for t < 0, all the particles acquired the 
average flow velocity jy in the x direction at t = and 
then the Lee-Edwards boundary condition Jl|,[l4J main- 
tained the simple shear flow. Steady sheared states were 
realized after transient viscoelastic behavior. 

As has been confirmed in the literature ]8[-|l0|], the 



dynamics in quiescent supercooled states is reasonably 
well described by the Rouse model, where the relax- 
ation time of the p th mode of a chain is expressed as 



(b 2 /[l2k B T sin 2 ( P Tr/2N)} (p 



!,••-, N-l) 



Here, the statistical segment length b is related to the 
variance of the end-to-end vector of a chain P — R^—Ri 
by 6 = [{P 2 )/{N-1)] 1 / 2 , which is 1.17, 1.18, 1.19 for T = 
1.0, 0.4, 0.2, respectively. We determined the monomer 
friction constant C from the relaxation (P(t) ■ -P(O)) = 
2N" 1 Efco.i,- cot 2 (7r(2^ + 1)/2N) exp(-t/r 2 f+i) of the 
end-to-end vector. The Rouse relaxation time tr(= T\ = 
(b 2 N 2 /3TT 2 k B T) then increases drastically with lowering 
T as t r = 250, 1800, and 6 x 10 4 for T = 1.0, 0.4, 0.2, 
respectively. We also calculated the a relaxation time 
T a from F q (T a ) = e _1 at q = 2n ||], where F q (t) = 

EJLi( ex P[*9 ■ i R j(t) - Rj{0)])/N is the van Hove self- 
correlation function for the displacement of a tagged par- 
ticle. Then we obtained the result t q = Q.QUQb 2 /k B T at 
any T, so we have tr = 2N 2 r a in our system. 

Now let us discuss the linear viscoelastic behavior in 
supercooled states. In Fig.l we show the stress relaxation 
function, 



G{t) = (a^(t)a^ v (0))/Vk B T, 



(1) 



where cj is the space integral of the xy component of 
the total stress tensor over the volume V = L 3 . At the 
lowest temperature T = 0.2, G(t) exhibits salient fea- 
tures of glassy polymer melts (jjfl]- Its initial value is 
of order 100 (in units of e/u 3 ) and is very large, and it 
relaxes to a value Go about 5 for t > 1. We then have a 
slow decay of the form, 



G(t)^G cM-(t/T s )% 



(2) 



where t s = 90 ~ r a and (3 — 0.5. The agreement to 
Eq.(2) is excellent for 1 < t < 10t s . This glassy be- 
havior arises from monomeric structural relaxation. For 
t > 50r s it approaches the Rouse stress relaxation func- 
tion, 



G R (t) = nk B TN 



N-l 

P =i 



cxp(-2i/r p ), 



(3) 



which is equal to T(N — 1)/^V for t < r Q and decays 
as TN^ 1 exp(— 2i/r^) for t ^> tr in the dimension- 
less units. Obviously, this final stage behavior arises 
from relaxation of large scale chain conformations. The 
crossover between these two regions occurs in a nar- 
row time time range in our case. Experimentally, how- 
ever, the intermediate region, which connects the glassy 
and polymeric (Rouse or reptation) relaxations, extends 
over a much wider time range (typically 4 decades [0), 
and G(t) there has been fitted to an algebraic form, 
G(t) = e- l G {t/T s )- a with a - 0.5 @,§. In addition, 
with increasing the molecular weight, a rubbery plateau 
has been observed to develop after the crossover before 



O 




FIG. 1. The stress relaxation function G(t) (thin-solid 
lines) at T = 0.2 in a supercooled state and T = 1 in a 
normal liquid state. It may be fitted to the stretched expo- 
nential form (dotted line) at relatively short times and tends 
to the Rouse relaxation function Gn(t) (bold-dashed lines) at 
long times. 

the terminal decay, whereas it is not apparently seen 
in our short chain system. In our case, the (zero- 
frequency) linear viscosity rj consists of a monomeric part 
Arj of order 10r s from the integration in the time region 
t < 10t s and the Rouse viscosity 77^ = L dtGu(t) = 
OMSTN' 1 ^ from t > tr. The ratio A?//^, is thus of 
order 1/(TJV) (~ 1 at T — 0.2), whereas we should have 
Arj ^C tjr for much larger N . 

In the Rouse model, the space integral of the polymer 
(entropic) stress cr?g is the sum of ksTbjabjp/b 2 over 
all the bonds in the system, where bj a are the Carte- 
sian components of the bond vectors bj = Rj+\ — Rj. 
We have confirmed that the relaxation function G c (t) — 
(crf y (i)a^ / (0))/V r fc s T nearly coincides with G R (t) = 
G(t) for t > O.lTfl, whereas it is about a half of Gn(t) for 
t ^5 T s- We note that the bond vectors have the nearly 
fixed length b = 0.96, and a bond orientation tensor Q a p 
may be defined as 



JV-l 



Qa^iN-ir^ibjab^/b 2 . 



(4) 



i=i 



Then we have Q a p oc (cr^o). If the electric polarizability 
tensor of a bead is uniaxial along the bond vector, the 
deviation of the dielectric tensor Ae Q/ a is proportional to 
Qq/3 — S a f3/3. In flow birefringence we have, 



Ae xy = C{a* y )/V, 



(5) 



where C is a constant. In supercooled states (c^ y )/V 
can be much smaller than the total shear stress a xy , for 
instance, in transient states or in oscillatory shear. 
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FIG. 2. Shear stress divided by shear rate j = 10 , 10 2 , 
10~ 3 , 10" 4 (thin-solid lines) vs t/r R (where tr = 6x 10 4 ) 
at T — 0.2. The curves follow the linear growth function 
(bold-solid line) for jt < 0.1, but afterwards depart from it. 
The linear growth function in the Rouse model is also plot- 
ted(bold-dashed line). 



Therefore, the usual stress-optical law Ae x 



Ca x 



which is valid far above T g: breaks down close to T g . 

In Fig. 2 we show the stress growth function a xy (t)/ A / 
after application of shear at t = at the lowest temper- 
ature T = 0.2 for various 7. The curves are the av- 
erages of data of eight independent runs. In the ini- 
tial stage, in which jt < 0.1, we can see the linear 
viscoelastic growth, cr xy {t)/j — L G(t')dt', whereas a 
nonlinear regime sets in for jt > 0.1, resulting in the 
non-Newtonian viscosity 77(7). As a guide, we also plot 
the linear growth function L Gft(t')dt' from the Rouse 
model, which is much smaller than the true linear growth 
for t -C Tr. The relevant physical processes are as fol- 
lows: For jt < 0.1 the overall chain conformations are 
affinely deformed, whereas for jt > 0.1 the structural re- 
arrangements among beads belonging to different chains 
become appreciably induced by shear (as in the case of 
supercooled simple fluids jy]). Experimentally, a stress 
overshoot (a rounded maximum of cr xy (t)) has been ob- 
served at jt — 0.05 — 0.1 for higher molecular weight 
melts close to T g fm. 

In Fig. 3 we display the steady state viscosity 77(7) at 
T = 0.2, 0.4, and 1. It exhibits marked shear-thinning 
behavior for jtr > 1 and becomes independent of T for 
very high shear rates. The horizontal arrows indicate the 
linear viscosity from the Rouse model tjr, and the ver- 
tical arrows indicate the points at which 7 = t^ . In 
particular, the curve of T = 0.2 may be fitted to r\ oc 7^" 
with v ~ 0.7 for 777? > 1. Similar shear-thinning has 
been reported in MD simulations of short chain systems 
in normal liquid states, but at much higher shear rates 
[16,Q~7|. In the case of supercooled simple fluid mixtures 
[11 1, shear-thinning becomes apparent for 7 > O.lr" 1 . 
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FIG. 3. The steady state viscosity vs 7 for T = 0.2, 0.4, 1. 
A line of slope —0.7 is also written as a view guide. 

In the present short chain system, significant shear- 
thinning occurs at much smaller shear of order t^ ~ 
T a 1 N~ 2 - At this early onset, the overall elongation of 
chains take place, as will become evident in Fig. 5 below. 
It is worth noting that the shear stress at 7 ~ t^ 1 is 
of order nksT/N (which would be the modulus in the 
rubbery plateau for longer chain systems). 

We finally examine anisotropy in the chain conforma- 
tions in shear at T = 0.2 and 7 = 10~ 4 . In Fig. 4 (a), we 
plot the xy cross section (z — 0) of the steady state bead 
distribution function, 



N 

g s (r) = N- 1 J2(S(Rj-RG-r)), 
3=1 



(6) 



where Ra{t) = N 1 Xm=i Rn{t) is the center of mass of 
a chain. In Fig. 4 (b), we plot the structure factor, 

N 

5(q) = iV- 2 ]T (exp(zq • (Rt - R,)), (7) 

in the q x % plane (q z = 0). It is proportional to the scat- 
tering intensity from labeled chains in shear E3|. We 
recognize that g s (r) and S(q) almost saturate into the 
forms in Fig. 4 for 7 > 10/tr. The figures indicate that 
our short chains take ellipsoidal shapes on the average 
once 7 > 10 /tr. The angle 9 between the ellipsoids and 
the y (shear gradient) direction is also written, which is 
calculated using Eq.(8) below. Let us define the tensor 

lap = EX=i(0R< " Rj)<*{Ri - Rj)l3)/N 2 - For small q 
with q z = 0, we have the expansion, 



S{q) = l- 1 -ai{q-e 1 ) 2 - l -al{q-e 2 f 



(8) 



where {ei,e2} and {a 2 ,a|} are the eigen unit vectors 
and values of the tensor I a p with a,0 = X,y. The two 
lengths a\ and ai are the shorter and longer radii of the 





FIG. 4. (a) Isointensity curves of <? s (r) in Eq.(6) in the 
xy plane (-3.75 < x,y < 3.75, z — 0). (b) Those of the 
incoherent structure factor S(q) in Eq.(7) in the q x q y plane 
(— 7r < q x ,q y < n, q z — 0). The values on the isolines are 
0.01 + 0.02n in (a) and 0.1 + 0.2n in (b) with n = 0, 1, • • • , 4 
from outer to inner. Here T = 0.2, j — 10~ 4 , and the flow 
is in the horizontal (x) direction. The 9 is the angle between 
the average chain shapes and the y axis. 
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FIG. 5. tan#, 1 — 01/02, and Q xy vs 7-77}. 



ellipses. In Fig. 5 we write tan = —e\ y /e\ x , the degree 
of elongation 1 — 01/02, and the xy component Q xy of 
the tensor in Eq.(4). For ttr > 10 we have 6 = 80°, 
ai/a 2 = 0.3, Qzy = 0.1. These quantities represent 
the average chain forms and bond orientation and are 
insensitive to T if plotted vs jtr. We stress that the 
shape changes of chains start to occur at 7 ~ r^ 1 while 
the monomeric structural relaxation is only slightly af- 
fected by shear. In fact, we found that the cage breakage 
time Tb(j) among neighboring beads belonging to differ- 



ent chains IO] does not change appreciably for 7 

at T = 0.2. This tendency should be more evident for 

longer chain systems. 

(i) In our simulations, the stress relaxation function 
obeys a stretched exponential decay and the Rouse re- 
laxation, although the chain length is too short and the 
temperature is too high to reproduce the glass-rubber 
transition region and the rubbery plateau, (ii) We have 



also found very early onset of the nonlinear regime of 
shear. Strong shear-thinning and anisotropic scattering 
can be predicted for 7 > r^ or for <j xy > nkgT/N in 
the Rouse case N < N e . In the entangled case N > N e 
the threshold shear rate and stress needed for the on- 
set of nonlinearity should be of order r~Jp and nksT/Ne, 
respectively, where r rep is the reptation time. Scatter- 
ing experiments from very weakly sheared melts near T g 
are promising, (hi) Although not presented here, hetero- 
geneities are much enhanced at low T in the cage break- 
age events among beads belonging to different chains p2| . 
They are characterized by the correlation length £ depen- 
dent on T and 7. Interesting crossover effects can then 
be expected at £ ~ bN 1 ' 2 for longer chain systems. 
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